function out = getAVarTheta11Step3Boot(theta11Step3,AVarTheta2Step2,AVarTheta12Step3,...
        setupStep1,setupStep3,outputStep3,optim,epsValue,tTestOn,rhoCrossAdj)
    
% Fixing the seeds
seedNum   = 1000;
rng(seedNum,'twister')

% required information for the bootstrap
numTheta11       = size(fieldnames(theta11Step3),1);
namesTheta12     = AVarTheta12Step3.names;
numTheta12       = size(namesTheta12,1);
posTheta12Step3  = zeros(1,numTheta12);
for i=1:numTheta12
    posTheta12Step3(1,i) = find(strcmp(AVarTheta12Step3.names,namesTheta12(i)));
end
theta11BootDraws     = zeros(setupStep3.numBootStep1And3,numTheta11);
tTestTheta11         = zeros(setupStep3.numBootStep1And3,numTheta11);
rhoCrossBoot         = zeros(setupStep3.numBootStep1And3,1);
yfitSample           = outputStep3.yHatTrans;
residualsTrans       = outputStep3.residuals';
[T,ny]               = size(residualsTrans);

% Recentering the residuals in the cross-section dimension 
for t=1:T
   select = ~isnan(residualsTrans(t,:));
   residualsTrans(t,select) = residualsTrans(t,select) - mean(residualsTrans(t,select),2);
end

% The bootstrapped indices used when resampling 
indexBoot        = zeros(setupStep3.numBootStep1And3,ny);
for b=1:setupStep3.numBootStep1And3
    indexBoot(b,:) = ceil(rand(ny,1)*ny);
end

% The AR(1) model to account for cross-correlation
estimRhoCrossOn = 0;
if setupStep3.wD > 1
    if isempty(rhoCrossAdj)
        rhoCross = BootstrapARone_CrossSection(residualsTrans);
        estimRhoCrossOn = 1;
    else
        rhoCross = rhoCrossAdj;
    end
else
    rhoCross = 0;
end
% EXTRA: for using the block bootstrap   
%indexBoot        = block_bootstrap([1:ny]',setupStep3.wD,setupStep3.numBootStep1And3);
%indexBoot        = squeeze(indexBoot(:,1,:))';
%dlmwrite('indexBoot.txt',indexBoot);    %For using the same index in Fortran

if setupStep3.wD > 1
    % we account for cross-sectional dependence. Missing values in residualTrans are then put to zero
    for t=1:T
       select = isnan(residualsTrans(t,:));
       residualsTrans(t,select) = 0;  
    end
end


%parfor (b=1:setupStep3.numBootStep1And3,2)  % For using multiprocessing with 2 workers
for b=1:setupStep3.numBootStep1And3         % For ordinary loop
    if mod(b,10) == 0
        disp(['Number of bootstraps so far = ',num2str(b)])
    end
    
    % Copy of the setup file
    setupStep3Boot = setupStep3;
    
    % Generate the bootstrap sample.
    eps = residualsTrans(:,indexBoot(b,:));  % standard resampling
    if setupStep3.wD > 1
        % we account for cross-sectional dependence
        for i=2:ny
            eps(:,i) = rhoCross*eps(:,i-1) + residualsTrans(:,indexBoot(b,i));
        end
    end
    setupStep3Boot.data = yfitSample + eps;
    
    % Update the values of Theta12 in the setup file.
    % We use the robust draws (if lambda = 0, then the robust and the homo draws are the same)
    for i=1:numTheta12
        name = namesTheta12(i,1);
        setupStep3Boot.calibrateTheta1.(name{1}) = AVarTheta12Step3.theta12DrawsRobust(posTheta12Step3(1,i),b);
    end
    
    % The optimization
    setupStep3Boot.optimizer = optim.optimizer;
    [~,theta11Step3Boot,outputStep3Boot] = step1SR(theta11Step3,setupStep3Boot,optim);
    theta11BootDraws(b,:) = struct2array(theta11Step3Boot);
    
    % Estimating the pooled AR(1) for cross-sectional dependence
    rhoCrossBoot(b,1) = BootstrapARone_CrossSection(outputStep3Boot.residuals');
    
    if tTestOn == 1
        % Standard errors of theta11
        AVarTheta11Step3Boot = getAVarTheta11Step3_usedInBootstrap(namesTheta12,theta11Step3Boot,AVarTheta2Step2,...
            setupStep1,setupStep3,epsValue,b);
        use_SE = sqrt(diag(AVarTheta11Step3Boot.VarHomo))';
        %use_SE = sqrt(diag(AVarTheta11Step3Boot.VarHomoPooled))';
        tTestTheta11(b,:)     = (theta11BootDraws(b,:) - struct2array(theta11Step3))./use_SE;
    end
end

%% Computing percentiles in the bootstrap sample
% The percentiles
prctilesValues     = [1.0 2.5 5.0 10 50 90 95 97.5 99];
theta11_pTileValues= zeros(numTheta11,size(prctilesValues,2)); 
for i=1:numTheta11
    theta11_pTileValues(i,:) = prctile(theta11BootDraws(:,i),prctilesValues);
end
% Storing the percentiles in a struct
for j=1:size(prctilesValues,2)
    label = num2str(prctilesValues(1,j));
    label = strrep(label,'.','_');
    theta11_pTile.(['pct',label]) = theta11_pTileValues(:,j);
end

%% output 
out.names                = fieldnames(theta11Step3);
out.theta11_se           = std(theta11BootDraws,0,1);
out.theta11_pTile        = theta11_pTile;
out.theta11_Draws        = theta11BootDraws;
out.nyObsBootstrap       = size(residualsTrans,2);
out.numBootStep1And3     = setupStep3.numBootStep1And3;
out.rhoCross             = rhoCross;
if estimRhoCrossOn == 1
    out.RhoCrossAdj      = 2*rhoCross - mean(rhoCrossBoot);    
    if abs(out.RhoCrossAdj) >= 1
        [~,sigma2,RvHatPooled] = BootstrapARone_CrossSection(residualsTrans);
        % We use: sigma2/(1-rhoCross^2) = RvHatPooled to determine rhoCross
        out.RhoCrossAdj = sqrt((1-sigma2/RvHatPooled));
        out.RhoCrossAdj_shrinkingOn = 1;
    end
end
if tTestOn == 1
    out.tTestTheta11 = tTestTheta11;
    % for tests which may not be symmetric
    tTestTheta11_pTileValues = zeros(numTheta11,size(prctilesValues,2)); 
    for i=1:numTheta11
        tTestTheta11_pTileValues(i,:) = prctile(tTestTheta11(:,i),prctilesValues);
    end
    % symmetric t-tests
    AbstTestTheta11_pTileValues = zeros(numTheta11,size(prctilesValues,2)); 
    for i=1:numTheta11
        AbstTestTheta11_pTileValues(i,:) = prctile(abs(tTestTheta11(:,i)),prctilesValues);
    end
    
    % Storing the percentiles nicely in a struct
    for j=1:size(prctilesValues,2)
        label = num2str(prctilesValues(1,j));
        label = strrep(label,'.','_');
        tTestTheta11_pTile.(['pct',label])    = tTestTheta11_pTileValues(:,j);
        AbstTestTheta11_pTile.(['pct',label]) = AbstTestTheta11_pTileValues(:,j);
    end
    out.tTestTheta11_pTile    = tTestTheta11_pTile;
    out.AbstTestTheta11_pTile = AbstTestTheta11_pTile;
end

% Turn-off multiprocessing
%delete(gcp)
end